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Abstract 

p/ ' The maintenance of large-scale differential rotation in stellar convective regions by ro- 

C/!5 ■ tationally influenced convective stresses also drives large-scale meridional flows by angular- 

i-^ ■ momentum conservation. This process is an example of "gyroscopic pumping", and has 

I ■ recently been studied in detail in the solar context. An important question concerns the 

i-i • extent to which these gyroscopically pumped meridional flows penetrate into nearby stably 

^ ■ stratified (radiative) regions, since they could potentially be an important source of non-local 

mixing. Here we present an extensive study of the gyroscopic pumping mechanism, using 
a combination of analytical calculations and numerical simulations both in Cartesian geom- 
etry and in spherical geometry. The various methods, when compared with one another, 
provide physical insight into the process itself, as well as increasingly sophisticated means of 
estimating the gyroscopic pumping rate. As an example of application, we investigate the 
i/^ ■ effects of this large-scale mixing process on the surface abundances of the light elements Li 

and Be for stars in the mass range 1.3-1.5Mq (so-called "Li-dip stars"). We find that gy- 
roscopic pumping is a very efficient mechanism for circulating material between the surface 
and the deep interior, so much in fact that it over-estimates Li and Be depletion by orders 
^ . of magnitude for stars on the hot side of the dip. However, when the diffusion of chemical 

species back into the surface convection zone is taken into account, a good fit with observed 
surface abundances of Li and Be as a function of stellar mass in the Hyades cluster can be 
found for reasonable choices of model parameters. 



1. Introduction 

1.1. Gyroscopic pumping 

Recently, Garaud & Acevedo-Arreguin (2009) (GAA09 hereafter) presented a prelimi- 
nary analysis of a "new" mechanism for rotational mixing in stellar interiors and conjectured 
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on its potential role in the depletion of Lithium in young Main Sequence stars. This mech- 
anism, called "gyroscopic pumping" , was originally studied in the context of Earth's atmo- 
spheric dynamics by Haynes et al. (1991) and later discussed in the astrophysical context by 
Gough & Mclntyre (1998) and Mclntyre (2007) who argued that it plays an important role 
in the dynamics of the solar interior. Loosely speaking, gyroscopic pumping occurs in any 
rotating fluid in the presence of additional stresses or forces which perpetually accelerate 
or decelerate the flows in the azimuthal direction. By conservation of angular momentum, 
the accelerated parcels of fluid move away from the rotation axis while the decelerated ones 
move toward the rotation axis, thus generating meridional fluid motion. 

Stars with outer convective regions often exhibit a significant amount of surface differen- 
tial rotation (e.g. Collier-Cameron 2007; Reiners 2007). This differential rotation is thought 
to be maintained by anisotropic and spatially varying Reynolds stresses (see Riidiger 1989, 
for example), which tend to continually accelerate the equatorial regions and decelerate the 
poles in a manner most remarkably observed in the solar convection zone (Schou et al. 1998). 
If one assumes that the star is close to dynamical equilibrium, its mean rotation rate f2^ lies 
in between the polar and equatorial rotation rates. When viewed in a frame rotating with 
angular velocity Q^, the differential rotation of the star's convective zone forms a large-scale 
azimuthal flow pattern, typically prograde in the equatorial region and retrograde in the 
polar regions. The aforementioned "gyroscopic pumping" can then be viewed in two equiv- 
alent ways. As described earlier the constant acceleration of the equatorial regions is a local 
source of angular-momentum to the fluid, which by angular-momentum conservation must 
move outward from the rotation axis. Similarly, fluid in the polar regions must move toward 
the rotation axis. Alternatively, one may simply note that the Coriolis force associated with 
the azimuthal flows described above pushes the fluid away from the rotation axis in prograde 
regions, and toward the rotation axis in retrograde regions. As illustrated in Figure [H the 
process naturally drives large-scale meridional flows throughout the outer convection region, 
with an upwelling near the equator, a poleward velocity near the surface and downwelling 
near the poles. The polar downwelling does not need to stop at the radiative-convective 
interface, and could in principle cause significant non-local mixing between the outer con- 
vection zone and the regions below. 

A similar gyroscopic pumping process is likely to drive fluid motion from within an inner 
convective zone as well. While the internal rotation profile of stars with convective cores has 
never been observed, one can readily expect some degree of differential rotation since the 
turbulent stresses associated with the rotationally constrained convective motions are likely 
to act in a similar fashion to those of the outer convective region. In that case again, there 
is no a-priori reason for the flows thus generated to stop at the interface with the overlying 
radiative zone, and one may wonder how much mixing they induce in the star. 
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Fig. 1. — Illustration of gyroscopic pumping in a star with an outer convection zone (shaded 
area). Slower rotating regions near the poles pump the fluid towards the rotation axis through 
the action of the Coriolis force (Fc) while the rapidly rotating region near the equator pump 
the fluid away from the rotation axis. A large-scale circulation is thus driven, poleward 
near the surface. The position of the return flow depends on the radiative zone physics (see 
GAA09 and gUfor detail). 

1.2. Mixing in radiative regions induced by gyroscopic pumping 



The first quantitative study of gyroscopic pumping in the context of stellar interiors 
was recently presented by GAA09. They focused on the solar case, i.e. a star with an outer 
convective region and an inner radiative region, and worked in the Boussinesq approxima- 
tion (arguing that the solar radiative zone does not span too many pressure and density 
scaleheights) . They showed that the fate of gyroscopically pumped flows - how much overall 
mixing they induce beyond the convective zone - depends equally on the thermal stratifica- 
tion and on the dynamical properties of the nearby radiative region. 

In accordance with the earlier results of Garaud & Brummell (2008), they found that 
in stratified, rotating stars in quasi-steady dynamical balance, the effect of stratification on 
large-scale meridional flows is principally controlled by the quantity 



a 



N 



(1) 



where Pr = u/k is the Prandtl number (where u is the local viscosity, and k is the local 
thermal diffusivity) , A^ is the local Brunt- Vaisala - or buoyancy - frequency and fi^ is the 
mean stellar rotation rate. When a is large, the effect of thermal stratification is strongly 
felt by the gyroscopically pumped meridional flows, which are exponentially damped away 
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from the radiative-convective interface on the lengthscale -R*/cr. Correspondingly, if a is 
small the flows can in principle penetrate much more deeply into the radiative interior, and 
cause significant large-scale mixing. 

It is important to note that a depends on fi^, so that the effect of stratification (in 
the sense defined above) is strongly reduced for rapidly rotating stars. Figure |2] shows an 
estimate of a for various stars in the mass range l.SM© — l.SM© at age 300Myr. In all 
cases, a remains well-below unity showing that large-scale mixing of the radiative zone by 
gyroscopic pumping could be significant for these stars. 
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Fig. 2. — Profile of a{r) for three different stellar masses at age 300Myr. The I.3M0, IAMq 
and I.5M0 stars are assumed to be rotating with the velocities 8 x 10~^rad/s, 1.3 x lO'^^rad/s, 
and 1.5 x 10~^rad/s respectively. The stellar models are computed using a standard stellar 
evolution code developed by Bodenheimer et al. (2007). 



Crucially, however, GAA09 showed that even in the very weakly stratified limit {a <^ 1), 
a second condition needs to be satisfied for the pumped flows to penetrate into the nearby 
radiative region. Indeed, this limit corresponds to the case where the system's dynamics 
are dominated by the balance between the perturbation to the pressure gradient and the 
Coriolis force. This well-known situation is called the "Taylor-Proudman" state, and (in the 
Boussinesq approximation) implies that all components of the velocity field are strongly con- 
strained to be constant along the rotation axis. GAA09 showed that the Taylor-Proudman 
constraint can prohibit flow generated within the outer convection zone from entering the 
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underlying radiative zone. Indeed, by mass conservation, any flows entering the radiative 
zone must somehow return to the convection zone. But such return flow would necessar- 
ily require breaking away from the Taylor-Proudman state. Hence, if there exists a region 
within the radiative zone where the Taylor-Proudman constraint is broken, then this region 
provides a channel through which the pumped flows can return, as illustrated in Figure [31 If 
such a region does not exist, the meridional flows instead return within the convection zone 
causing negligible mixing in the radiative zone. 

Since the Taylor-Proudman constraint is broken whenever there exist additional stresses 
of amplitude comparable with the Coriolis force, various mechanisms can be invoked. Of 
particular interest are Maxwell stresses in the presence of small- or large-scale magnetic fields 
(see Gough & Mclntyre 1998, for example), and turbulent stresses within another convective 
region. It is the latter we are mostly interested in studying in this work, namely the case of 
stars with a convective envelope and a convective core. 



Q. 




Fig. 3. — Illustration of the expected flow pattern for rapidly rotating stars with two convec- 
tive zones. Gyroscopically pumped flows from one convection zone can enter the radiative 
zone and return within the second convective zone where the Taylor-Proudman constraint 
does not apply (see main text for detail). Note the role of the tangent cylinder in such a 
model in delimiting mixed regions from non-mixed regions. 



This paper is organized as follows. We flrst present a fairly exhaustive study of gyro- 
scopic pumping for rapidly rotating stars (in the limit where a < 1). In ^we lay out the 
general formulation of the problem and the assumptions made. Note that the basic model 
used in this paper is largely inspired from the work of GAA09 but extends it to the case 
of multiple convective regions and uses a more generally applicable formalism for the fluid 
dynamics, namely the anelastic approximation. In ^ we flrst solve the problem analytically 
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in a Cartesian coordinate system, using a much simplified stellar model. This exercise pro- 
vides insight onto how the meridional flow velocities induced by gyroscopic pumping scale 
with the forcing mechanism (in this case the differential rotation of the inner and/or outer 
convection zone), with the stratification within the radiative region, and with the system's 
geometry (i.e. the respective widths of the inner and outer convective regions). In ^ we 
then apply the same model in a two-dimensional spherical geometry, using a more realistic 
stellar model as the background state. By comparing Cartesian-model results with equiva- 
lent spherical-model results, we deduce a very simple rule to go from one to the other. This 
rule is particularly useful since one-dimensional Cartesian model numerical solutions can be 
obtained in a tiny fraction of the time necessary to integrate the two-dimensional spherical 
case, and can also be pushed to true stellar parameter values (which cannot be done in two 
dimensions) . 

Next, we present a simple application of this theory to estimate the surface abundances 
of lithium (Li) and beryllium (Be) in young Main Sequence stars as a function of their age and 
mass. More precisely, we are interested in young stars in the mass range of 1.3Mq — 1.5Mq, 
which have the well-known property of being significantly depleted in Li and Be in their 
surface layers as compared with slightly more and slightly less massive stars (Boesgaard & 
Tripico 1986; see Boesgaard 2005 and Anthony-Twarog et al. 2009 for reviews). These so- 
called "Li-dip" stars are unique in the sense that they have two significant convection zones 
(one inner and one outer) which, as described earlier, could promote large-scale mixing 
between the surface and the interior by gyroscopic pumping. In §3 we show that the Li and 
Be depletion rates as induced by gyroscopic pumping for Li-dip stars can be very significant. 
For reasonable model assumptions, the cool (low-mass) side of the dip is readily explained by 
gyroscopic pumping, while depletion fractions on the hot (high-mass) side are much larger 
than observed. When the diffusion of chemical species back into the outer convection zone 
(by overshooting motions for example) is taken into account, good agreement between the 
model and the data is achieved on both sides of the dip. Finally, we conclude in ^ by 
summarizing our main results and discussing future prospects. 



2. The model 

In this section we briefly derive the model equations used throughout this paper. We 
consider a star of radius R^, rotating with a mean angular velocity Q^,. In all that follows, 
we assume that the system is in a quasi-steady state, non-magnetic, and axially symmetric. 
While these assumptions are probably over-simplistic, they help us narrow down gyroscopic 
pumping to its essence. 
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The star considered can have up to two convective regions: a convective core, extending 
from the center of the star to a first radiative-convective interface located at the radius 
r = Tin; and a convective envelope, located between the outer radiative-convective interface 
at r = Tout and the surface. Note that the stellar surface and both interfaces are assumed to 
be perfectly spherical. 

GAA09 studied large-scale meridional flows within the solar interior using the Boussi- 
nesq approximation. This approximation treats the thermodynamical quantities as the sum 
of a weakly varying background plus small perturbations. One may argue in favor of its 
use in the solar radiative zone, which does not span too many density scaleheights. Here, 
however, we aim to model stars with very thin outer convective regions, in which case the 
radiative zone extends nearly all the way to the stellar photosphere. We must therefore 
switch to using the more general anelastic approximation instead, which allows for a more 
strongly varying background stratiflcation. 

The anelastic approximation implicitly assumes that all velocities are small compared 
with the local sound speed, and that thermodynamical perturbations are small compared 
with the equivalent background quantities. We thus deflne p{r), p{r), T{r) and s{r) as 
the spherically symmetric background pressure, density, temperature and entropy proflles, 
and the equivalent p, p, T and s as two-dimensional perturbations to these quantities. If the 
equation of state is assumed to be that of a perfect gas (which is an acceptable approximation 
for the purpose of this work), then the perturbations are related by 

p p T 

neglecting for simplicity the dependence on the chemical gradients. 

The momentum, thermal energy and mass conservation equations describing the dynam- 
ics of the interior flows and thermodynamical perturbations, in the anelastic approximation, 
are: 

Vv P I 
u ■ Vu + 2fl^ X u = — — + -g + -V ■ n , 

P _ P P 
pcpfu . Vs = V ■ (krVf) - V ■ Ft , 

V ■ (pu) = , (3) 

where u is the velocity fleld expressed in a frame rotating with angular velocity Q^,, g = 
(0, 0, —g) is gravity, 11 is the viscous stress tensor, Cp is the speciflc heat at constant pressure, 
/ct = pCpK, is the thermal conductivity and flnally Ft is the turbulent heat flux (in the 
convective regions). Note that we have neglected any distortion of the star caused by the 



centrifugal force, as well as perturbations to the gravitational field. These assumptions 
suppress the well-known (c/o 6a/ Eddington- Sweet flows (see Spiegel & Zahn 1992). 

Following GAA09, we model the inertial term of the momentum equation in the convec- 
tive regions by the sum of a turbulent viscosity plus a linear drag term driving the system 
toward a differentially rotating profile: — (u — 'Ucze0)/T, where MczG,/, is the assumed/observed 
azimuthal velocity profile in the convective zones and where r is the local convective turnover 
timescale (which varies with depth). This drag term is introduced to "mimic" the effect of 
turbulent convection on driving differential rotation. It is only significant in the convection 
zones, and drops to zero in the radiative region. The very slow flow velocities expected 
in the radiative zone justify neglecting the various nonlinear terms in the momentum and 
heat equation there. We replace the turbulent heat advection term by a turbulent diffu- 
sivity, which is assumed to be very large in the convection zones and rapidly tends to zero 
otherwise. 

The resulting model equations, which we use throughout this work unless otherwise 
specified, are therefore: 

2f2, X U = ^ + ^g + _V ■ (n + Hturbj 



p p p ^ 

^ u-e, = V-((fcT + A;turb)VT) 

9 



1 = P^Z 
p p T 

V • (pu) = , (4) 

where Ilturb is similar to the microscopic stress tensor, but using a turbulent viscosity in- 
stead. Note that we have rewritten the background heat advection term to emphasize the 
dependence on the buoyancy frequency N (see Spiegel & Zahn 1992, for example). 



3. A simplified Cartesian model 

Much can be learned about the dynamics of stellar interiors by first studying a simplified 
problem in Cartesian geometry (the "planar star" approximation, see Garaud & Brummell 
2008 and GAA09 for example). Equations in this geometry can usually be solved analytically 
to gain insight into the physical processes at play. They often reveal important scaling laws 
governing the solutions, and finally, are rarely more than an "order one" geometrical factor 
away from more realistic solutions in spherical geometry (in fact we prove this in ^. Our 
primary goals in this section are therefore not quantitative. Rather, we aim to determine. 
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qualitatively, how deeply the meridional flow velocities generated by gyroscopic pumping 
penetrate into the radiative zone, and characterize how their amplitude scales with the 
system parameters. 

Since the Cartesian model solutions obtained are knowingly off by a factor of order 
unity anyway, we further simplify the equations, in this section only, with the following 
substitutions: 

- V ■ (n + Hturb) ^ (/> + //turb)V2u , 
P 

V ■ {{kT + A;turn) Vr) ^ {It + A;turb) V'T . (5) 

It can be shown (through numerical integrations) that neither of these substitutions affect 
the scalings of the solutions in the case where both convection zones are presenlo. Finally, 
and following Spiegel & Zahn (1992) we neglect in this section, for analytical simplicity, the 
pressure perturbations in the linearized equation of state so that: 

p T 

p T ^ ^ 



3.1. Model setup and non-dimensional equations 

As in GAA09, we consider a Cartesian coordinate system {x, y, z) with z aligned with 
both gravity and with the rotation axis. The x-direction represents the azimuthal direction, 
while the y-direction is equivalent to minus the co-latitude. Distances are normalized to 
the stellar radius i?*, so that the stellar interior is in the interval z G [0,1], the inner 
convection zone spans the interval [0, z-^^^ and the outer convection zone spans [-Zout, !]• Figure 
m illustrates the geometry of the Cartesian system. We assume that the star is axially 
symmetric, i.e. independent of x, and periodic in y on the interval [0, tt] (representing the 
two "poles") with equatorial symmetry (i.e. symmetric about y = vr/2). 

The background stellar model is chosen to be very simple, so that analytical solutions 
of the problem can easily be found. We take 

• p(z) = peC"^/^" and T{z) = T^e-'/^^, 



^The substitution of the stress tensor, however, modifies the nature of the viscous boundary layers (the 
well-known Ekman layers, see for example Kundu 1990), which are relevant in the case of stars with a single 
convection zone only. 
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z=l 



z=h„ 



z=hi, 



z=0 



y=Jl/2 



Fig. 4. — Cartesian model geometry and intended correspondence with the spherical case. 
The shaded area marks the convective regions, where forcing is applied. The y = and 
y = 7r/2 lines mark the "poles" and the "equator". The system is assumed to be periodic 
with period tt in the j/— direction. 

• u, R and g are constant. 



Note that fits to actual stellar models show that the non-dimensional density and temperature 
scaleheights typically satisfy the inequalities Dp <^ Dt < 1. 

The global velocity field in the rotating frame is u = {u, v, w) and flow velocities are 
normalized to -R*f2^. In this framework, the unit timescale is fi^^. Density and temperature 
perturbations are normalized by pc and eT^ respectively, where e = R-^Q^/g is the ratio of 
the centrifugal force to gravity. Pressure perturbations are normalized to pcR^^I- The set 
of equations dl]) then simplifies to the non-dimensional system: 






U- Un 



r 



f^? 



V • fe-^/^^u) = , 



where we have introduced a series of standard parameters, namely the Ekman number 

V 



E. 



' Rl^. 



and the equivalently defined -E'l/.turb, as well as 



E. 



K 



Rin, 



(7) 



(8) 



(9) 
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which is actually an inverse Peclet number, and the equivalently defined i^K^turb- Note that 
the microscopic diffusivities normally vary with depth within a star, but are assumed here to 
be constant for simplicity. Deep within the interiors of Hyades-age stars in the Li-dip mass 
range. 



E,, 



10^^^- 10 
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10" 



(10) 



In order to fully specify the model, the functions N{z), t{z), Ucziy, z) and the turbulent 
diffusivity profiles must be selected. The buoyancy frequency is taken to be 



N{z) 



2 



tanh 






+ tanh 



z- z-„ 



Ai 



;ii) 



so that iV(2;) — )> in both convection zones, and iV(z) ~ N^^ in the radiative zone. The 
lengthscales Ain and Aout Kiay be thought of as the respective thicknesses of the "overshoot" 
regions located near each of the two convection zones (and in this section are taken to be 
equal to one another, for simplicity). For t(^z\ we take 



/ \-l ^*-in 

Az) = -IT 



tanh 



A, 



A 



out 



tanh 






;i2) 



where Ain and Aout are assumed to be constant and equal to the inverse of the (non- 
dimensional) convective turnover time in the relevant convection zone. 

Both convective zones may be differentially rotating. For mathematical simplicity again, 
we assume in this Cartesian model that their latitudinal dependence is similar. We therefore 
select the following functional form for Ucz{y,z): 



Ucziy,z) 



UrAz e 



iky 



(13) 



where k = 2 to guarantee equatorial symmetry. The radial profile Ucz{z) is then chosen to 
be 



u°!^Hz) 



UcAz 



1 + tanh 



Z Zq^i 

^out 



Kz[z) 



1 + tanh 



Aw 



(14) 



The functions u°^^{z) and u'^{z) can a priori depend on depth (see GAA09 for example). By 
analogy, the non-dimensional turbulent diffusivity profiles are constructed as 



E: 



out / 
!^,turb \ 



EuturhyZ) 



and similarly for Ei^^turh{z). 



tanh 



Z Zq^i 

^out 



E, 



U,tUTh \ 



tanh 



Z]ri Z 



Ai 



(15) 
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Projecting the model equations into Cartesian coordinates, using invariance in the 
X— direction and seeking periodic solutions in the form of q{y, z) = q{z)e^^'^ for each of 
the dependent variables yields: 



-2v = {E^ + ^i.,turb) [u^^ - k'^u) - 



U — Ur 



T 



2u = -ikpe''^" + {E^ + E^,turb) [vzz - k^v) - - , 



w 

T 



= -p.e^/^" + Te^/^^ + {E, + E,,turb) {w,, - k^w) - 

_e-./i?T^ = [E^ + E,,turb) (t,z - k^T^ 

ikve-'/^' + {e-'/^'w), = . (16) 

where the subscript z denotes a derivative with respect to z. 

Finally, we need to specify an adequate set of boundary conditions for the system. The 
two boundaries at 2 = and z = 1 are assumed to be impermeable {w = 0), stress- free 
{uz = Vz = 0), and the temperature perturbations are assumed to be zero. Note that as long 
as the system boundaries are located in a convective region, the actual choice of boundary 
conditions has little influence on the result. 



3.2. Solution and interpretation of the model 

The set of equations (fT6l) and associated boundary conditions (see above) can be solved 
analytically when the overshoot regions are very thin compared with the depths of the 
respective convective zones. The complete derivation of the solution is fairly straightforward 
although algebraically cumbersome. It is detailed in Appendix A: exact solutions are derived 
in each of the three regions [0, z^^], [zi^, -Zout] and [zout, 1], and matched to one another across 
the radiative-convective interfaces at z-m and Zout respectively. In what follows, we discuss 
the most important outcome of this analysis, namely the prediction of the gyroscopically 
pumped mass flux mixing the radiative zone. 



3.2.1. General behavior of the model solutions. 

As found by GAA09, in the limit where a < 1 the meridional flows generated in the 
convective regions can penetrate deeply into the nearby radiative zone. A very simple way 
of seeing this is to note that in this limit, the x— component of the momentum equation in 
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f lT6|l reduces to t) = 0{E,,) in the radiative zone which then imphes, by mass conservation, 
that 

(e-^/^^w), = 0{E,) . (17) 

Hence the non-dimensional vertical mass flux mixing the radiative zone, pw, is constant 
along the rotation axis: 

pw = W,, , (18) 

and spans the entire region, extending from one convection zone to the other as drawn in 
Figure [3] for example. 

The details of the calculation of the pumped mass flux Wrz? even for this simplifled 
stellar model, are fairly complicated and are presented in Appendix A. In the limit where 
the depths of both convection zones d\-a = -Zin and liout = 1 — -Zout are small compared with 
the stellar radius (which is true for most stars in the Li dip), we show that 

2 P 



Wr: 



Q2 Ek k^ \ Aindin ~^ Aoutdout 



+ 



where P is the "pumping term" 



P(^outXz (2:out) - piZin)u'^^iZin) + d, 



d{pu 



out\ 



out 



dz 



+ 4 



d{pul 



dz 



(19) 



(20) 



and where G is the following fairly obscure factor: 



Li — L (2;out ~ Z: 



+ kL^ 
where 



111/ 
2L 2L 



~dout/L 



'^out Wout 

2L 2L 






Zin/L 



sinh(/c) 

{Zout-Zin)/L 

sinh(A;) 



-'out 

IT 



3-(2out-2in)/i 



1 + 



sinh(fc) 

''' '' d- 



L ) kL 



sinh(/c) 



1 
kl 



(21) 



L-i = D-^ - D^^ . (22) 

Note that G is found to be negative for all reasonable parameter values, so that the denom- 
inator of flT9|) never vanishes. Our analytical solution is easily verifled by comparison with 
numerical solutions of the full set of equations (ITB|) . as shown in §3.2.31 However, let us flrst 
attempt to understand the meaning of ( IT9ll on physical grounds. 



3.2.2. Interpretation of the dependence of the pumped mass flux on physical parameters 



This expression for PVrz calculated in §3.2.11 can be interpreted more easily in two dif- 
ferent asymptotic limits. 
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The unstratified limit. When A^rz — > (the unstratified hmit), Wrz takes the simpler 
foriio: 

l^-l = It^a^^ITX^ ■ (23) 

r\j in I out 

A-inrfin Aoutdout 

We see that PKz is roughly of the order of the pumping term, times a factor which depends 
only on the respective properties (depth, convective turnover time) of the convective zones. 
Our main conclusion is that the mass flux into the radiative zone appears to go to zerqj if 
one of the convection zones vanishes, either by becoming vanishingly thin [din or dout -^ 0), 
or if the associated convective stresses become negligible (Ajn or Aout — ^ 0). 

This important result can in fact be easily understood in the light of the work of GAA09 
described in §1.21 Indeed, in the unstratified case, radiative zone flows must satisfy the 
Taylor-Proudman constraint. Any flows generated by gyroscopic pumping in one convection 
zone can only enter the radiative zone if there is a return path at the other end. If the second 
convective zone vanishes, this return path is no longer available. The flows instead return 
within the existing convection zone, and do not mix the radiative region significantly. The 
stratified case with a < 1 is very similar. 

Based on these very simple considerations, we can therefore expect that the overall 
mixing rate in the star resulting from gyroscopic pumping must reach a maximum for a 
given stellar mass between IMq (no or negligible convective core) and I.SMq (no or negligible 
convective envelope). This simple idea motivated our study of the Li dip (see ^, although, 
as we shall show, the real problem is much more subtle. 



The stratified case. The unstratified limit discussed above is of course artificial. In real 
stars A^rz 7^ and one must instead compare the two terms in the denominator of ( !T9|l to 



^Note that this expression for Wrz can be derived directly, and much more easily, by considering an 
unstratified system in the first place (ignoring the buoyancy term in the momentum equation, taking p and 
T constant, and ignoring the thermal energy equation). 

■^In practice, if one of the convection zones vanishes entirely (din = or dout = 0) then the analysis 
presented in Appendix A is no longer valid. It can be shown instead (analytically and numerically) that the 
meridional flow amplitudes do not entirely drop to zero but instead drop to the level of Ekman (viscous) flows 
and depend sensitively on the boundary conditions. Meanwhile, if Ajn = or Aout = 0, but the turbulent 
stresses £^"turb ^^'^ ^""urb ^^'^ non-zero, then mixing of the radiative zone by the pumped flows can still be 
effective. Indeed, gyroscopic pumping from one of the two convective zones is still effective, and the second 
provides the return pathway for the flows. This effect is not expressed in (J19p . since our analytical derivation 
ignores for simplicity the effect of the turbulent diffusion term compared with the relaxation term. 
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one another. In the hmit where the first term is much larger than the second then 

\W I - ^^^ (24) 

The flow velocities pumped into the radiative zone are now found to follow a local Eddington- 
Sweet scaling law, which is not surprising since we are looking at quasi-steady flows in a 
stratified fluid driven by rotational forcing. Such solutions were already found by Spiegel & 
Zahn (1992) for example in the case of the Sun. Naturally, local Eddington-Sweet flows are 
much slower than the flows pumped through each individual convective zone, although they 
could still provide significant sources of mixing in fairly rapidly rotating stars. 

The effect discussed in the unstratified case, namely the complete suppression of W^z in 
the limit where one of the convection zone disappears, still occurs but only when 

^ = ^^ « 1 (^±3l + 4 + AL \ .25) 

E^ ni E^ A;2 V AinC^in Aoutt^out / ■ 

Since N^^/Ef^Ql is typically very large because E,^ is very small, this limit is only relevant 
for extremely thin or weak convective regions. 



Based on these considerations, we can now re-interpret (1191) in the following way: W^z 
has two strict upper limits: a first upper limit, which arises from mechanical constraints (by 
the driving force and the Taylor-Proudman constraint), and a second upper limit which arises 
from the thermal stratification of the system (which can strongly suppress radial flows). The 
actual value of l^rz mixing the radiative zone is the smaller of the two, thus defining two 
different regimes, the "unstratified" regime (where Wj-^, is mechanically constrained) and the 
"weakly stratified" regime (where W,.z is thermally constrained). 



3.2.3. Comparison with full solutions of the governing equations 

In this section we verify the analytical solution for W^z expressed in ( [T9l) by comparison 
with numerical solutions of the same equations. Numerical solutions are obtained by inte- 
grating the two-point boundary value problem ( IT6|) with associated boundary conditions, 
background state and forcing as specified in §3.2.11 

A direct comparison of the analytical and numerical solutions for W^z is shown in Figure 
[5l for a wide range of simulations. For all the calculations presented, the microscopic diffusion 
parameters E^ = 10^^ and the Prandtl number Pr = E^/E,^ = 10"^ are fixed. The turbulent 
diffusivities are set to for simplicity. For the forcing by the differential rotation, as expressed 
in flT^ . we take u°^^{z) = ^^."(-z) = mq and Ain = Aout = 10. Note that the problem is linear 
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so the value of uo is irrelevant. The overshoot depths are taken to be Ain = Aout 



10- 



The background density and temperature scaleheights are selected to be Dp = 0.075 and 
Dt = 0.275, so that L = 0.1. Note that these scaleheights are a good approximation to the 
true density and temperature scaleheights in the radiative zones of stars in the 1.3Mq — 1.5Mq 
range. The non-dimensional buoyancy frequency of the radiative zone N^^/Vt^, is varied 
for each set of simulations, leading to a ranging from 10~^ to 10. Finally, each symbol 
corresponds to a particular geometry of the system, with different possible pairs of values 
(^in, ^^out) as shown. 

For each simulation the quantity pw is measured from the numerical solution at z = 0.5. 
Note that for a < 1, pw is found to be constant across much of the radiative zone (as 
expected), so exactly where pw is measured does not matter much. For cr > 1, pw is no longer 
constant, so the specific height z = 0.5 is selected for consistency across simulations. We 




Fig. 5. — Comparison between the numerical and analytical results for W^z = pw. For each 
selected value of a ranging from 10~^ to 10, two sets of simulations are shown with different 
convection zone depths: the case of thinnest convection zones {z^-^ = 0.05, Zo^t = 0.95) 
is shown as diamond symbols, while (^in = 0.15, z^^t = 0.85) is shown as squares. The 
corresponding analytical solutions are shown as dotted and solid lines respectively. For a < 1 
the numerical and analytical results fit very well for all cases, showing that the approximate 
formula for Wrz (see equation (iT9l) ) applies. For a > 1 the formula doesn't do so well, as 
expected. 



find that the analytical solution flT9|) . shown in the dotted and solid lines, fits the numerical 
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results very well as long as a < 1. The mismatcli for a > 1 was expected because we can 
no longer assume that pw is constant across the radiative zone (see GAA09 for example), 
while this is key to the derivation of f ll9p . Figure E] also shows the two regimes discussed in 
§3.2.3[ In the "unstratified limit", for a — )■ 0, we observe a plateau in W-^^ which sets the 
largest achievable value for flows in the interior. As the degree of stratification increases, the 
flows driven into the radiative zone are slowed down by the local thermal stratification, a 
phenomenon which becomes more pronounced as N^.^ increases (equivalently, as a increases). 
The flow velocities in this "weakly stratified" regime scale as o"~^, which is equivalent to the 
aforementioned local Eddington-Sweet scaling since E^j and Pr are constant. Note that in 
this regime, the flow pattern still spans the entire region with constant pw. 



3.3. Summary 

To summarize this section, we have used a simple toy model to study how gyroscopic 
pumping by convective zone stresses induces large-scale meridional fluid motions. We found 
that the pumping indeed drives significant large-scale flows within the convection zone(s) 
(see GAA09). Furthermore, in some circumstances, a fraction of the pumped mass flux may 
enter the adjacent radiative zone, and induce a circulation of material with the following 
properties: 

• Convection zone flows penetrating into the nearby radiative zone are exponentially 
damped on a lengthscale -R^/u where a is given by equation ([T]). If a < 1 (which is 
the case for most rapid rotators) the flows can potentially mix the entire star (Garaud 
& Brummell 2008, GAA09). 

• In the limit a < 1 the pumped mass flux within the radiative zone pw is mechanically 
constrained to be constant along the rotation axis. 

• The fraction of the gyroscopically pumped mass flux which does enter the radiative 
zone (by contrast with the mass flux which returns immediately within the driving 
convection zone) is capped by the lowest of two constraints: a mechanical constraint, 
which crucially depends on the presence of another source of stresses somewhere else 
within the system to enable flows to return to their point of origin (see previous sections 
for detail, as well as GAA09), and a thermal constraint, which limits the flow velocities 
to local Eddington-Sweet velocities. Note that this second constraint was not discussed 
by GAA09, but turns out to be the most relevant one for most stars. 
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The efficiency of gyroscopic pumping on mixing stellar interiors thus depends on many 
factors, including the background stellar structure, the thermal diffusivity, the stellar rotation 
rate and the nature of the convective zone stresses. However, thanks to the analytical formula 
fin?]) , we now have a reasonably clear picture of precisely how all of these factors influence 
the amplitude of the gyroscopically pumped mass flux within a star. 



4. From Cartesian to Spherical models 

The next step of this investigation is to move to more realistic numerical simulations 
of the problem in a spherical geometry. We have two goals in this endeavour. The first is 
to understand the effect of the spherical geometry on the solutions. In particular, we are 
interested in the dichotomy between the regions located respectively within and outside of 
the cylinder tangent to the convective core and aligned with the rotation axis (see Figure |3]). 
The second goal is much more quantitative, and is to extract (if possible) simple laws relating 
the calculated mass flux Wj-^ in the Cartesian case to its equivalent in the spherical case. 
Since Cartesian solutions are much easier to calculate than spherical geometry solutions, a 
simple rule to go from one set of solutions to the other could prove particularly useful later. 



4.1. Model description 

The spherical model used is very similar to the model presented in GAA09. The salient 
points are repeated here for completeness. 

We consider a spherical coordinate system (r, 6, (p) where ^ = denotes the rotation axis 
and 6 = 7t/2 marks the equator. The governing equations are the original ones (jl]) derived 
in §21 Since we are principally interested in studying the effects of the spherical geometry on 
the model predictions for the large-scale flow amplitudes, we consider "hypothetical stars" 
instead of real stellar models. This largely facilitates the comparison between the various 
simulation outputs, and enables us to focus on how the model depends on specific control 
parameters. 

In all cases, the "star" used is a solar- type star (i?^ = Rq, M^ = Mq), and the back- 
ground thermodynamical quantities such as density, pressure, and temperature (p(r), p{r) 
and T{r) respectively) are extracted from Model S of Christensen-Dalsgaard et al. (1996). 
However, to model cases with various size convection zones, we create artificial proffies of 
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the buoyancy frequency by using the expression 



N^(r) =Nlsm 



TT r — n, 



'"out ^ir 



if r G [rin, rout] 



N^(r) = -10 



-11 



otherwise, 



(26) 



where rin is the radius of the lower radiative-convective interface, and rout is the radius of 
the upper radiative-convective interface. The maximum value of N{r) within the radiative 
zone, A^rz; is one of our input parameters, while that of the convection zones is merely chosen 
to be arbitrarily low and has little effect on the outcome of the simulation. In what follows, 
we define the global parameter a^ as 



(J* 



n. 



(27) 



Note that the Prandtl number is assumed to be constant (see below). 



In all simulations, the star is assumed to be rotating at the mean angular velocity 
fi^ = 3 X 10~^rad/s. In all cases, the applied forcing by differential rotation is 



Ur 
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with 



^cz(6') = f^cq (l — 0.2 COS^ O) 



02 = 0.01 and figq = fi^ ( 1 — 7^ 



3a2 
IE 
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(28) 
(29) 

(30) 



to ensure that the total applied angular momentum to the system is zero. Note that for 
simplicity, the same differential rotation is chosen in both convection zones, and is measured 
by the parameter 02. This value is chosen to be fairly small, since the observed differential 
rotation of rapidly rotating stars is typically quite small. In practice, since we are studying 
a linear problem, the solutions scale linearly with 02. Finally, the expression for the non- 
dimensional quantity T{r) is the same as that given in ( !T2l) with Ajn = Aout = 10 in both 
inner and outer convection zones. 

The total diffusivities (i.e. the sum of the microscopic and turbulent components) are 
assumed to have the following profiles: 



K-{r) + K,tuvh{r) 



K^e^l^'^ + 



^turb 



2 + tanh ( — - — ^ | + tanh 



A. 



Ai 



y{r) + t'turb(r) = Pr (/t(r) + fi:turb(^)) 



(31) 
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The selected exponential profile for the microscopic part of h{r) is not too dissimilar from 
that of a real star in the 1.3-1.5Mq range provided D^^/R^, = 0.14. We define the inverse 
Peclet number E^c = i^c/Rl^*- The value of Ei^^ will be varied in the various simulations, 
and decreased as much as possible to reach the asymptotic stellar regime. We take the 
Prandtl number to be constant and equal to Pr= 10~^. The selected value of Pr is chosen 
to be smaller than one to respect stellar conditions, but larger than the actual stellar values 
(which are of order of 10~^ typically) to ease the numerical computations. The "convective" 
value Kturb is fixed so that the non-dimensional -E^.turb is equal to 10~^. The overshoot layer 
depths Ain and Aout are taken to be 0.01. While these choices are fairly arbitrary, they are 
reasonable given our qualitative goals. 

The computational domain is a spherical shell with the outer boundary located at 
0.95-R^ and the inner boundary at O.Oli?^,. The outer boundary is chosen to be well-below 
the stellar surface to avoid numerical complications related to the very rapidly changing 
background in the region r > 0.95/2^. The inner boundary is chosen to be well-within the 
inner convection zone, but excludes the origin to avoid coordinate singularities. The upper 
and lower boundaries are assumed to be impermeable and stress-free, with T = 0. 

The numerical method of solution is based on the expansion of the governing equations 
onto the spherical coordinate system, followed by their projection onto Chebishev polyno- 
mials Tn{cos6), and finally, solution of the resulting ODE system in r using a Newton- 
Raphson-Kantorovich algorithm. The typical solutions shown have 3000 meshpoints and 70 
Fourier modes. For more detail on the numerical algorithm, see Garaud (2001) and Garaud 
& Garaud (2008). 



4.2. Typical solution 

Solutions have been computed for a wide range of values of the parameters di^ = Tin and 
c^out = R* — ''^out (the respective depths of the inner and outer convection zones), S^c and a^. 
For low enough diffusivities the overall structure of the radiative zone flows converges to a 
pattern which only depends on di^, dont and a-^ (and with an amplitude which scales as the 
calculated Wrz). Figure [6] shows a representative example of the kind of 2D flow structure 
found within the star, for E^^^ = 10"'^ and cr^ = 0.1. An artificial case with fairly thick 
convective zones {r^^ = 0.2i?^ and rout = 0.8-R^) was chosen to make it easier to visualise the 
results. 
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Fig. 6. — Numerical solution for rout = 0.8i?^, r^ 
quadrants represent the stellar interior, with the radiative-convective interfaces denoted as 
dotted circles. The left-side panel shows selected streamlines of the flow, as logarithmically- 
spaced contours of the stream-function. Solid lines represent clockwise flow, while dotted 
lines represent counter-clockwise flow. The right-side panel shows the normalized differential 
rotation profile. 



Discussion of the meridional flow structure. Representative streamlines are shown 
in the left-side panel of Figure [6] and reveal the structure of the meridional flows driven 
by gyroscopic pumping. In the outer convection zone, we observe a dominant cell in mid 
and high latitudes, poleward near the surface and equatorward near rout- In addition, a 
small counter-cell of fairly surprising shape is observed in the equatorial region. The inner 
convection zone also has a dominant cell of the same vorticity (with poleward flows in the 
outer layers, and a deep equatorward return flow), and a small counter-cell located just above 
it along the polar axis. The structure of the two dominant cells in the respective convective 
zones are easily understood from angular momentum balance. The smaller counter-cell 
above the inner convection zone is required to match the flows downwelling from the outer 
convection zone to the inner core flows. 

The bulk of the mass flux generated in a given convection zone returns within, or close 
to the edge of that same convection zone. However, weak flows from the polar region of the 
outer convective zone, and from the equatorial region of the inner convection zone escape and 
mix the cylinder tangent to the convective core. Note how, by contrast, the region outside 
of the tangent cylinder is mostly quiescent. 

The variation of the vertical flow amplitude with latitude for the simulation of Figure 
E] can be seen more clearly in Figure [71 it illustrates the downward pumping of the flows in 



-22- 



the polar regions, while mixing in the equatorial region is mostly negligible. The variation of 
the amplitude of the flows with depth is well-explained by the variation in the background 
density: flows pumped downward along the polar axis roughly satisfy Ur = W^^/p (where Ur 
is the radial velocity) where Wj-^ is constant. This constraint is more easily seen in Figure [S] 
which shows the profile of pUr at 88° latitude for the simulation shown in Figure [61 Within 
the radiative zone we see that pUr is roughly constant as expected. Figure [8] also illustrates 
how only a small fraction of the pumped mass flux enters the radiative zone, while most of 
it remains within the generating convection zone. 







20 
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Fig. 7. — Variation of Ur with latitude for four different non-dimensional radii r/R^ for the 
simulation shown in Figure El The latitudinal position of the tangent cylinder at each radius 
is marked by a small diamond. This figure illustrates the polar downwelling, the return flow 
just outside the tangent cylinder and the mostly negligible mixing at low latitudes. Note 
that the calculated value of Ur should not be taken at face value: it needs to be rescaled 
appropriately in order to estimate mixing velocities in real stars (see §3.2.ip . See ^for more 
quantitative estimates in Li-dip stars. 



Discussion of the azimuthal flow structure. As seen in Figure Ej both convection 
zones exhibit a rotation profile close to the imposed profile (i.e. with Q/Q^q ranging from 
1-02 to 1 between the pole and the equator, and with Q nearly constant with radius), as 
expected. Meanwhile, the radiative zone exhibits a similar level of differential rotation, with 
a striking shear layer near the tangent cylinder. This shear is presumably caused by the 
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Fig. 8. — Variation of p\ur\ with depth near the polar axis (at a latitude of 88°) in the 
simulation of Figured A logarithmic scale is used to visualize simultaneously the convection 
zone flows and the radiative zone flows. The solid line denotes upwelling flows, while the 
dashed line denotes downwelling flows. Note the radiative region r/i?^ G [0.2, 0.8] where pUr 
is constant as expected from the analytical solution. 

deposition of negative angular momentum by the meridional flows as they carry fluid away 
from the rotation axis and begin to rise up in the radiative interior again. Other numerical 
simulations (not shown here) show that this feature is stronger when flows are stronger, i.e. 
when the system is more weakly thermally stratified or more rapidly rotating (a* smaller). 
An interesting consequence of this shear layer, however, is the possibility that it may become 
unstable to Rayleigh instabilities in the vicinity of the tangent cylinder, where for extreme 
cases the specific angular momentum may locally begin to decrease with distance from the 
polar axis. 



4.3. Comparison between Cartesian model predictions and spherical model 

solutions 



The results presented in the previous section show that, for stars with cr^ < 1, the 
downward pumped mass flux into the radiative zone W^z is indeed constant with depth along 
the rotation axis, as found in the Cartesian model analysis of ^ We can now compare more 
quantitatively the results of these spherical simulations with equivalent Cartesian solutions 
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for VFrz, first in order to verify the scalings derived in flT9|) and then to see if there exists a 
simple relationship between the Cartesian model velocities and the spherical model velocities. 

We first run a series of two-dimensional spherical simulations, based on the model pre- 
sented in §4.H with the following parameters varied: 

• We consider three different geometries: {vi^/R^ = 0.05, rout/R-k = 0.9) {vi^/Ri, = 
0.1, Tout/R* = 0.8), and (rjn/i?^ = 0.2, Vout/R* = 0.8). The purpose is to explore the 
effect of varying the convection zone sizes on the predicted velocities. 

• For the case with rin/-R* = 0.1, Vont/R* = 0.8 we consider two different values of cr^, (by 
changing A^rz): cTj, = 0.1 and a^ = 0.03. For the other two geometries a^, is fixed to be 
0.1. 

• Finally, we consider a range of inverse Peclet numbers from 10~^ down to the lowest 
achievable value, E^^^ = 10^^°. 

In all cases, we measure the mass fiux W^^^'^'^ = p\ur\ at r/_R^ = 0.5, at a latitude of 85°. 
Note that this choice is fairly arbitrary: p\ur\ does not change with depth nor with latitude 
much as long as the point selected lies well-within the tangent cylinder. 

In order to compare the value obtained in the spherical case with Cartesian model sim- 
ulations, we integrate the equivalent equations and boundary conditions, now expressed in 
a Cartesian coordinate systeno using exactly the same geometries {z-m = nn/i?*, Zout = 
Vout/R*), background profiles, diffusivities and convective turnover timescale as in the spher- 
ical cases described above (e.g. equation (fT2|) . and equation (13T!) with r replaced by z). We 
construct the forcing velocity Uc7,{z)e^^'^ based on the differential rotation profile ilcz(^, ^) in 
the following way: we take k = 2, and set 



Ucz{z) = — 



2 + tanh ( — - — ^ 1 + tanh 



An„t / V A 



(32) 



in both inner and outer convection zones. We then measure W^^"^^ = pw at the same height 
(2 = 0.5). 

As shown in FigureOwe find that there is an excellent agreement of the Cartesian model 
calculations with the spherical model calculations, for all sets of simulations at low enough 



"^Note that these equations are different from the ones used in ^and are indeed the Cartesian expression 
of (j4]). Specifically, they differ from those of [3] by using the correct viscous stress tensor, the correct heat 
flux, and the full linearized equation of state. 
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values of the diffusivities, provided we divide the Cartesian model results by a factor of 2tt: 

lycart 

The discrepancy for higher values of the diffusivities appears to arise when the diffusive layer 
thicknesses become of the order of the modelled structures (i.e. the thickness of the outer 
convection zone or the width of the tangent cylinder). 

This result implies that the overall scalings derived are indeed correct, but furthermore 
that it is possible to use the simplified Cartesian model to get very precise estimates of the 
mass flux into the radiative zone within the tangent cylinder. The good fit between the 
two sets of simulations can presumably be attributed to the fact that as long as Tin <^ i?^,, 
the tangent cylinder is quite thin and curvature effects should indeed be negligible. The 
multiplicative factor of 1/2tt is not obvious a priori (hence the need for this exercise), but is 
not particularly surprising either. We attribute it to the fact that in the cylindrical case, the 
rotation axis is "infinitely far" away from the region where the flows are calculated, whereas 
it plays an important role in the calculation in the spherical case. 



4.4. Summary 

To summarize this section, we have seen that the gyroscopic pumping mechanism still 
works (as expected!) in a spherical geometry, and that the dynamics within the cylinder 
tangent to the inner convective core (for stars with two convective zones, see Figure [3]) 
are very similar, qualitatively and quantitatively, to the dynamics studied in §3] In par- 
ticular, we found that the pumped mass flux within this tangent cylinder, as measured in 
full two-dimensional spherical geometry calculations, is equal to the equivalent pumped flux 
calculated in the Cartesian case, but divided by a factor of 27r. 

This important result provides an interesting and practical mean of getting precise 
estimates for the rate of mixing induced by gyroscopic pumping, which can in principle 
be used in stellar evolution models. In the following section, we provide an example of 
application of this idea, to the Li-dip problem. 
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Fig. 9. — Comparison between p\ur\ extracted from the spherical model calculations (sym- 
bols) as described in the main text, and the equivalent Cartesian model calculation (lines) 
where W^^'^^ is divided by a factor of 27r. The solid lines correspond to the case where 
Tin = 0.1-Ry, and rout = 0.8, for cr^ = 0.1 and c^ = 0.03 for the lower and upper line respec- 
tively. The dashed line corresponds to rjn = 0.2/?^ and rout = 0.8/?^, while the dotted line 
corresponds to rin = 0.1i?+ and rout = 0.9-R+. 

5. Application to the Li dip problem 
5.1. Introduction 



An outstanding problem in stellar astrophysics concerns the measured abundance of the 
rare light element Lithium in the atmospheres of F-type main-sequence stars (for reviews 
see Boesgaard 2005 and Anthony-Twarog et al. 2009). Lithium burns by nuclear reactions 
at temperatures of 2.5 xlO^ K or above in these stars, but the depletion is not evident 
at the surface unless there is a mixing mechanism to bring the Li down to layers at that 
temperature at some time during the evolution of the star. In the mass range of interest, 
between 1.1-1.6 M©, the pre-main- sequence convection zone does not extend down to high 
enough temperatures to result in appreciable Li depletion at the surface, and in fact most 
main sequence Pop I stars in this mass range have Li abundances {Nu ~ 2 x 10~^ that of 
hydrogen by number) characteristic of those in the youngest stars, indicating little, if any 
depletion. 
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However, as first discovered in the Hyades cluster (Boesgaard & Tripico 1986) there is 
a narrow range in effective temperature 6400 K < T^s < 6900 K (spectral types F6-F0) 
where a sharp dip in the Li abundance is observed, with a minimum value of Nu < 10^^^ 
at Tcff ~ 6650 K. The mass range within the dip is 1.3-1.5 M©. The depth of the surface 
convection zone decreases rapidly as Tgfj increases across the dip. The same temperature 
range also corresponds to a rapid change in spectroscopic rotational velocities (Boesgaard 
1987; Wolff & Simon 1997), with the stars around T^s = 7000 K rotating with vsini up to 
150 km/s and those at Teg = 6400 K with only 20 km/s. A similar dip is also observed in 
the Praesepe cluster (Soderblom et al. 1993a) with an age similar to that of the Hyades 
(600-700 Myr). In the much younger Pleiades cluster (100 Myr) the dip at about T^s = 6700 
K is marginal (Soderblom et al. 1993b) or not present (Boesgaard 2005). In the even 
younger cluster a Per (50 Myr) the dip is also not yet evident (Balachandran et al. 1996). 
Nevertheless the data suggest that at least some Li depletion occurs relatively early, before an 
age of 200 Myr (Anthony-Twarog et al. 2009). In the Hyades, a similar dip, but not as deep, 
is observed for the light element beryllium (Boesgaard & King 2002). In the older cluster 
IC 4651 (1-2 Gyr) both the Li dip and the (less deep) Be dip are observed (Smiljanic et al. 
2010). For further details on the properties of the dip in various clusters see Pinsonneault 
(1997) and Anthony-Twarog et al. (2009). 

Main-sequence surface convection zones in this mass range do not extend deep enough 
to mix Li and Be down to their respective burning radii, so the challenge is to find another 
mixing process that operates only in this particular range of spectral types. As summarized 
by Pinsonneault (1997) and Anthony-Twarog et al (2009), the various proposed mixing 
mechanisms to explain the dip can be divided roughly into three types: mass loss, diffusion, 
or slow mixing as a consequence of rotation or waves. Schramm et al. (1990) suggest that the 
temperature range of the lithium dip also corresponds to that of the pulsational instability 
strip where it intersects the main sequence, so that a slow mass loss rate, induced by low- 
amplitude pulsations, could simply remove the lithium remaining in the surface layers. 

Michaud (1986; see also Richer & Michaud 1993) explained the dip by diffusion and 
gravitational settling of Li atoms out the bottom of the convection zone. In their model the 
cool side of the dip arises from the increasing effectiveness of diffusion once the convection 
zone becomes thin, and the hot side is explained by radiative upward acceleration which 
counteracts the diffusion once the star becomes hot enough. To obtain good agreement 
with observations, a small amount of mass loss is also required in the theory. However, 
diffusion models tend to deplete Li and Be at about the same rate, and are not consistent 
with observations. 

Some form of rotationally-induced mixing seems to be the most promising effect to ex- 
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plain the observations (Pinsonneault 1997), in particular the Li/Be ratio in the dip (Dehyan- 
nis & Pinsonneault 1997). This process can circulate Li out of the surface convection zone 
down to layers where it can be destroyed, but no entirely satisfactory model has yet been 
found. Such mixing can be induced by gravity waves generated by the surface convection 
zone (Garcia Lopez & Spruit 1991, Talon & Charbonnel 2003) or by meridional circulation 
or secular shear instabilities (Deliyannis & Pinsonneault 1997). In the rotational mixing 
models the cool side of the dip is explained by the increase in rotational velocity as Tg^ in- 
creases, possibly combined with the gravity wave model. The hot side is much more difficult 
to explain; Talon & Charbonnel (1998) propose a model involving wind-driven meridional 
circulation and turbulent transport induced by differential rotation, based on earlier work 
by Zahn (1992) and Talon & Zahn (1997). This type of model was shown to be consistent 
with both Li and Be observations around the gap in IC 4651 (Smiljanic et al. 2010). 



5.2. Li and Be depletion by gyroscopic pumping 

In this section, we are primarily interested in determining the effect of gyroscopic pump- 
ing on Li and Be depletion as a stand-alone mechanism (i.e. in the absence of any of the effects 
described in §5.ip . We consider stars in the mass range of the Li dip, namely 1.3Mq — 1.6Mq. 
We use the results of [0] to estimate the depletion rate of Li and Be in the surface layers of 
these stars, induced by gyroscopic pumping, as follows. First, we find numerical estimates 
for the pumped mass fiux out of the convective envelope and fiowing into the deep interior 
within the cylinder tangent to the inner core (see Figure [3] and also Figure [TT!) . In order 
to do this, we solve the set of equations (jl]) expressed in a Cartesian geometry using a real 
stellar background model, extract the desired value of the pumped mass fiux W^^^^, and then 
use the rule (1331) to estimate the equivalent pumped mass fiux in the more realistic case of a 
spherical star. Using simple geometrical arguments, we then construct and solve evolution 
equations for the surface Li and Be abundances, which can be compared with observations. 



5.2.1. Background model 

We use the code developed by Bodenheimer et al. (2007) to construct a sequence of 
reference background models in the range 1.3M0 — I.6M0, evolved from the ZAMS up to 
300 Myr which is about half the age of the Hyades cluster. For reference, this code solves the 
standard equations of stellar structure and evolution, and uses a gray model atmosphere as 
an outer boundary condition. We assume a solar initial composition. The code is calibrated 
to match the Sun's observed properties at 4.57 Gyr. 
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Table 1 summarizes various properties of these modelled stars at age 300 Myr: the 
stellar radius R-^, the respective depths of the inner and outer convection zones d^^ and 
dont, an estimate of the convective velocities in the bulk of each convection zone, v^^^^ and 
v°Q*y and finally the radii ru and vbs below which the local stellar temperature exceeds the 
Li-burning and Be-burning temperatures of 2.5 x 10^ K and 3.5 x 10^ K respectively. 

These reference stars are assumed to rotate with an angular velocity i7^ derived from 
the results of Wolff & Simon (1997), who provide estimates for the mean vsini of stars 
in various mass ranges and ages (see their Table 4). We first note that in the mass range 
considered, the mean rotational velocities do not change much between the Pleiades age 
and the Hyades age. We interpolate the observations of these two clusters to an age of 
approximately 300 Myr. We also note that at the Pleiades age their rather high 1.4 — I.5M0 
measurement only has 4 data points; we discard it. We then interpolate their remaining 
results to our selected stellar masses to get: < vsini >obs (1-3M0) = 50km/s, < t>sini >obs 
(I.35M0) = 61km/s, < vsini >obs (1-4M0) = 79km/s, < vsini >obs (I.45M0) = 96km/s, 
< f sini >obs (1-5M0) = 98km/s, < vsini >obs (I.55M0) = lOOkm/s and < vsini >obs 
(I.6M0) = lOOkm/s. To convert the mean vsini measurements into rotational velocities, 
we note that the mean value of sin 2 over all possible measurements is (very roughly) 2/tt. 
In that case, we take 

2R^ ^ ' 

The resulting f2^ values are listed in Table 1. Note that because of the change in the stellar 
radius across the selected mass range, f2^ does not vary too much. 

We also extract from the stellar models the radial density, temperature, opacity and 
buoyancy frequency profiles within the stars, which are used as the background state for 
the Cartesian calculation of the gyroscopically pumped mass flux. From these quantities, we 
continue constructing the model as follows. The microscopic viscosity and thermal diffusivity 
profiles ^{r) and R{r) are calculated using the formula given by Gough (2007) (see also 
Garaud & Garaud 2008). The turbulent diffusivities in the convective regions (see equation 
f fTSj) ) are derived from the model convective velocities as: 

-^™turb(^) = -^"turb(^) = 020 ' ^"^^^ 

in the convective core, and similarly for the outer convection zone. The overshoot depths 
Ain and Aout niay be different near the inner and outer convective regions, and are taken 
to be equal to 10% of the local pressure scaleheight at the respective radiative-convective 
interfaces. 

The convection zones are both assumed to be rotating differentially with the profile 
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given in equation fl29|) . The parameter 02, which to a good approximation is equal to the 
difference between the equator and the polar rotation rate, normalized by fi^, is assumed for 
simplicity to be the same in both convection zones, and is a free parameter in the problem 
(recall that the predicted velocities scale linearly with 02). The convective velocities are 
used to derive the quantities Ain and Aout which characterize the relaxation timescale to this 
assumed differential rotation profile (see equation (TT2l) ). We take 

and similarly for Aout- Note that this quantity is the most difficult to relate to real stellar 
parameters, since it refers to our very simplified parametrization of the effects of turbulence. 
Here we have used R^, as a typical lengthscale instead of the depth of the convective region 
in calculating the relaxation timescale. The reasoning behind this choice is that angular 
momentum has to be transported all the way from the pole to the equator (i.e. a typical 
lengthscale i?*) for the large-scale differential rotation profile to be established. 

Finally, it is crucial to note that while some of our choices (of Q^,, of Ajn and Aout) are 
arguably arbitrary, it so happens that they do not influence the resulting depletion rates 
much. The reasons for this will be explained in detail in §5.2.31 



5.2.2. Calculated mass flux and depletion timescale 

In order to evaluate the mass flux induced by gyroscopic pumping we now integrate the 
two-point boundary value problem (jlj), expanded in Cartesian geometry, with the realistic 
background stellar model described in the previous section. The advantage of the Cartesian 
calculation is that it can indeed be performed with true stellar values of the diffusivities, 
using 100,000 meshpoints. This would not be possible with spherical-geometry calculations. 
We extract from the numerical solutions the quantity W^^'^^ by measuring the value of pw at 
z = 0.5, although, as discussed in §21 the exact height does not matter since pw is constant 
within the radiative zone. We first convert this value back into dimensional form, and then 
to the spherical case using (!33|) to obtain an estimate for W^^^'^'^, the local mass flux flowing 
in the tangent cylinder. The results are shown in Figure [TOl and discussed in §5.2.31 

Based on the numerical results from §4.H and best illustrated in Figure El we now 
approximate the Li circulation pattern in the manner depicted in Figure [TT] (see also Figure 
|3|). Li-rich fluid is pumped down from the outer convective zone, through the cylindrical 
region Cu (delimited in radius by the cylinder tangent to the convective core and in the 
vertical direction by ru and rout respectively), down to the Li-burning interior. By mass 
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Fig. 10. — Calculated mass flux W^^^'^^ in the model simulations as a function of stellar mass. 
Note that V^^p'^'^'' scales linearly with a2, hence the rescaling on the y— axis. To estimate a 
radial velocity at radius r, divide W^^^^^ by the density at that radius, p(r). 



conservation, Li-free material is pumped up from the deep interior, through Cu, to the 
outer convective zone. As a result, Cu and the outer convection zone are both progressively 
depleted in Li with time. 

Note that the timescale for material to flow vertically from the Li-burning radius ru 
to the base of the outer convection zone rout (or vice-versa) can easily be calculated by 
integrating the inverse of the radial velocity W^^^'^^/p between ru and rout- The results, 
using W^l^*^^ estimated above and the true stellar density proflles, are shown in Figure [T21 
together with the equivalent timescale for the transport of Be from the Be-burning radius 
r-Be to the outer convection zone. In all stars considered, the transport timescales are of 
the order of 30-100 Myr and 100-300 Myr for Li and Be respectively when the differential 
rotation rate a2 is one percent. Note that the transport timescale is not equal to the depletion 
timescale, although the two are related (see Appendix B for detail). Figure ITU also shows 
the similarly calculated timescale required for material to move from the convective core to 
the convective envelope. It is interesting and reassuring to note that for reasonable values of 
02, this timescale is of the order of tens of Gyr, so that one would not expect to see a strong 
modiflcation of the surface abundances of He or CNO nuclei over the lifetime of these stars. 
This is an important self-consistency check of the model. 
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Fig. 11. — Schematic of the Li (or Be) mass flux in the various regions of a Li-dip star's 
interior. The surface and the two radiative-convective interfaces are shown as solid circles. 
The dashed hne represents the Li-burning (or alternatively the Be-burning) radius - only one 
of the two is shown on this Figure for clarity. Fluid is pumped from the surface downward 
along the rotation axis, and returns to the outer convection zone in a thin layer close to the 
edge of the tangent cylinder (vertical dashed line). By mass conservation, both fluxes are the 



same and equal to W^^^*^ 



Tirt 



The mass abundance of Li (and Be) in the outer convection 



zone is X^'^, that in the cylinder CLi,Be is X'^^^ while X = below the burning radius. In the 
region below the outer convection zone but above the burning radius, the abundance {Xq) 
is essentially primordial. Diffusion (squiggly arrows) brings Li and Be back into the outer 
convection zone. 



We now evaluate the rate of change of the mass fraction of Li in the outer convective 
zone, denoted as X£j . By definition X£? = Ml?/M°"*, the ratio of the total mass of Li in the 
outer convection zone to the total mass of the outer convection zone. Since material flows 
through Cu on its way up and down, we also need to deflne the mass fraction of Li in Cu-, 
Xl^ = M^f /M{Cu), i-e. the total mass of Li contained in Cu divided by the mass of that 
region. Note that 
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r'out rsin ^(rin/r) 
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(37) 



and can easily be integrated numerically. 
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Fig. 12. — Timescales for the circulation of material from the base of the outer convective 
zone down to the Li-burning radius (solid line), the Be-burning radius (dashed line) and the 
inner convective core (dotted line). Note that since W^^^'^'^ scales linearly with 02, so do these 
timescales, hence the rescaling of the y— axis. For reasonable values of 02 (of the order of 
a percent), the Li- and Be-circulation timescales are smaller than the age of the star, while 
the circulation from the surface down to the inner core is of tens of Gyr. 



Based on Figure [HI we see that 
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since the total mass flux both up and down, within the tangent cylinder, is equal to 
Tirf^W^^^'^'^ . Dividing the first equation by M°^^ and the second by M{Cu), we finally get 
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which implicitly defines two timescales. 
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The timescale Tp^ is the characteristic timescale over which the material within the outer 
convection zone is recirculated by gyroscopic pumping, while T^p^^^ is the timescale over 



which the material within the cylinder Cu is recirculated. If t^'^^^ ^ '^if pump then X^? ~ X^f 

'"pump is of the order of or greater than Tl^pu^p 



at all times. Meanwhile, if Tp^ is of the order of or greater than TLf , then X£; can 



differ from X^f significantly. 

The set of equations (139!) can easily be solved analytically given the initial condition 
X£?(t = 0) = Xl^ (t = 0) = Xufi, the initial Li abundance in the star. The calculation 
is detailed in Appendix B. A very similar calculation can be done to estimate the surface 
Be abundance Xg^ as a function of time. The only difference is that one should replace r^ 
with TBe, the Be-burning radius, which defines the slightly larger cylindrical region Cse, and 
equivalently the quantity Xg^ . The timescale Tp^^j^p remains the same, while we define 

Be,pump ^^2^yf/_spher " I ) 

Our theoretical results, now expressed as Li and Be depletion fractions, are shown in 
Figure [T3l assuming values of 02 = 0.002, 0.005, and 0.01. They are compared with observed 
Li and Be abundances in the Hyades as reported by Boesgaard (2005). It is quite clear 
that contrary to our naive expectation of [|3l the gyroscopic pumping of Li and Be out of 
the outer convection zone does not decrease for the higher-mass stars despite the fact that 
their outer convection zone becomes smaller and smaller. This is also clear from Figure [T0| 
where W^^^'^^ continues to increase with stellar mass instead of being quenched as initially 
expected. In the next sections, we discuss why this is the case, and how to reconcile the 
model with observations. 



5.2.3. Discussion 

The most important conclusion from this analysis is that gyroscopic pumping alone may 
be able to explain the cool side of the dip for reasonable values of the differential rotation 
rate (one merely needs to adjust 02), but always vastly over-estimates the depletion rate for 
stars on the hot side of the dip. Since the gyroscopic pumping mechanism is quite generic, 
and arises from simple first-principles of angular momentum conservation our study then 
raises the question of how one might suppress the effects of pumping in the hot side of the 
dip. In order to address the problem, it is important first to understand the cause of the 
sharp decrease in the predicted depletion timescale for high-mass stars in this model. 

First, note that most of the pumping in the stars considered comes from the inner core. 
Indeed, the density in the outer convection zone is so low that the mass flux pu°"* generated 
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Fig. 13. — Relative depletion of Li and Be in the Hyades cluster (assumed to be 700Myr 
old) as a function of stellar mass around the Li dip. Note that AA^(Li) = log(X£;/XLio) 
and similarly for Be. The thin lines connecting the symbols show data interpolated from 
Boesgaard (2005). The solid line and squares correspond to the Li data, while the dashed 
lines and triangle correspond to the Be data. The model predictions are shown in thick solid 
and dashed lines for the Li and Be cases respectively, from top to bottom in each case for 
02 = 0.002, 0.005 and 0.01 (the larger 02, the larger the depletion). 



in the outer convection zone is negligible compared with the mass flux generated from the 
inner convection zone (see equation (fT9|) ). As a result, the characteristics of the inner core 
dominate the forcing of the flows, while the outer convection zone merely plays the role of 
providing a pathway for the flows to return to the interior. This is easily verified numerically 
by setting (artificially) Aout = to suppress the forcing in the outer convection zone. Since 
turbulent viscosity in the outer convection zone is nevertheless still present, the return path 
still exists and as seen in Figure [TU the resulting depletion rates are hardly changed. 



One may then wonder what the main factor controlling the variation in the depletion 
rate as a function of stellar mass actually is. A few immediate possibilities come to mind. 
The convective velocities in the inner core of these stars, as well as the size of the inner core 
(see Table 1), both increase with M^, implying that the total mass flux pumped by the inner 
convection zone is larger for stars on the hot side of the dip. In addition, the increase in the 
rotation rates of the stars as M^ increases implies that the allowed flow speeds through the 
radiative zone are larger on the hot side of the dip (see equation (I24p ) . All of these effects 
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explain the trend seen in Figure [TOl which clearly shows that the pumped mass flux W^l^*^^ 
increases with M^ across the dip. 

However, this is not sufficient to explain the vast increase in the depletion rates observed 
in Figure [13] as M^ increases. In Figure [H] we show different artificial models to illustrate 
this statement. Model Ml is the aforementioned case where Aout is set to zero. Model M2 is 
created holding fl^, constant and equal to lO'^^rad/s across all stars in the model, to suppress 
the effect of increased rotation rate across the dip. Model M3 is created holding Ajn constant 
and equal to 10 across all stars in the model, to suppress the effect of increased convective 
velocities across the dip. Model M4 is created holding Tin constant and equal to 0.05i?* 
across all stars in the model, to suppress the effect of increased core size across the dip. In 
all cases, all other quantities are the same as the original model of Figure [T31 As we can 
see, none of these changes affect the predicted depletion rates much. This incidentally also 
shows that the exact details of the model are not particularly important, and that another, 
much more fundamental effect controls the overall depletion rate. Finally, we show (as the 
two dotted lines) a variant of model M2 in which Q^, is held constant but this time equal to 
2 X 10~^rad/s. This value is closer to the typical angular velocity (derived from vsini) of 
stars for which Li has actually been observed, which is significantly smaller than the median 
rotation rates for stars of the same mass (see Boesgaard 1987). Very similar depletion rates 
to those of models M1-M4 can be recovered provided the overall differential rotation is chosen 
to be larger (02 = 0.08). This shows that there is some degeneracy in the model parameters, 
which is not entirely surprising. 

The dominant effect in the model depletion trend is in fact found to be the decrease in 
the mass of the regions which need to be depleted in Li or Be as M^ increases. This mass 
has two contributions: the mass of the outer convective zone M^^^, plus the mass contained 
in the cylinders Cu and Cse respectively. To illustrate the effect in question, we create two 
additional artificial models. In model M5 the mass of the outer convective zone is artificially 
held constant (and equal to 1O~^M0). In model M6 the mass in the cylinders is set to 
(assuming that only the Li and Be fractions in the convective zone need to be recirculated). 
The results are also shown in Figure [HI The depletion rates in each case are now strikingly 
different from those of models Ml-4. In the case of M5, the depletion fraction is now much 
more constant across all stars. This strongly suggests that the depletion timescale in Figure 
[13] varies with stellar mass more because the total mass which needs to be depleted varies 
than because the pumped mass flux W^^^^"^ varies. Model M6 illustrates a fundamental 
property of these rapidly rotating Li-dip stars undergoing gyroscopic pumping. Since the 
pumped mass flux is independent of depth below the convection zone, unless the mass of 
Be to be depleted is different from the mass of Li to be depleted, the predicted depletion 
timescales will be exactly the same for the two species. This is precisely the case illustrated 



-37- 



o.c 


^"W--.. ..„„= 


-'.0 


; \ % - 


_ ■• r 


\ % : 


-2.0 


\ \ ^ 


-2.5 


- 


^ 


-3.0 




- 



1.3 



1.4 



M5 



=.:;!= M1-M4 



Ml 
M4 
M3 
M2 



.6 



MVM^ 



Fig. 14. — Relative depletion of Li (solid lines) and Be (dashed lines) in the Hyades cluster 
as a function of stellar mass in the Li dip, as predicted in various artificial test-models. AA^ 
is defined as in Fig. [131 In all cases, 02 = 0.002. In model Ml, Aout is set to zero for all 
stars. In model M2, fi^ = 10^'^rad/s for all stars. The dotted lines show the variant of model 
M2 for which Q^ = 2 x lO^^rad/s and 02 = 0.08, for Li and Be respectively. In model M3, 



Ai- 



10 for all stars. In model M4, r-,. 



0.05R^ for all stars. In model M5, the mass of 



the outer convective zone is set 1O~^M0 for all stars. In model M6, the contribution of the 
cylinder (see main text) to Mu and M^e is set to 0. All other quantities are as in Figure [T3l 
Note that in model M6, the two lines for Li and Be are on top of each other, as expected 
from the analytical model. 



by model M6 which ignores the contribution of the cylinders Cu and Cbb, so the predicted 
depletion fractions of Li and Be are exactly the same. Hence, the mixed inner cylinders are 
crucial to the difference in the depletion timescales of Li and Be. Moreover, their radii and 
heights uniquely determine the predicted ratio of the Li and Be depletion fraction. As seen 
in Figure [131 the original model described in §5.2.21 appears to predict the correct ratio, even 
though the absolute depletion fractions are too large for the hot side of the dip. 

The remaining question is why the pumping is not suppressed for stars of masses higher 
than I.5M0, as originally expected from the gradual disappearance of the outer convective 
zone. In fact, it turns out that c?out and v°^^^ remain significant in these stars, even though 
the total mass of the convection zone becomes negligible. Since the main role of the outer 
convection zone is to provide a pathway to recirculate the small mass flux pumped by the 



-38- 



inner core, there is no notable suppression of the pumping in this mass range contrary to 
our original idea described in ||3l 



5.3. The role of Li (and Be) diffusion 

While the model as it stands appears to reproduce the depletion trend on the cool 
side of the dip with reasonable assumptions for the differential rotation of the inner core 
(02 = 0(0.01)), it vastly overestimates the depletion fraction on the hot side of the dip. Our 
efforts now shift to the problem of reconciling the model with observations, or in other words 
on reducing the depletion rate on the hot side of the dip. 

The answer, as it happens, is quite simple, and lies in the balance between advection of 
Li-rich (and Be-rich) material out of the outer convection zone by gyroscopic pumping, and 
diffusion of these elements back into it from the radiative zone below which, as illustrated in 
Figure fTTl still has primordial Li and Be abundances. To take this new effect into account, 
equation f l38|) must be modified as 

M- = -nrlw:t'^Xn + 7rrlW:t'^X^' - A7Trl,p{r,^,)DVXu , (42) 

where the new third term is the mass flux of Li diffused back into the outer convection 
zone. The diffusion coefficient D is presumably the sum of a microscopic and a turbulent 
contribution from convective overshoot. The microscopic contributioiu is evaluated using 
the prescription of Gough (2007), while the turbulent contribution presumably decays ex- 
ponentially with depth below the convection zone on a typical overshoot depth lengthscale, 
yielding 

D{r) = Dmicro(r) + Aurb(r) = ^micro(r) + V°^^^dout exp ( ) for r < T^nt ■ (43) 

The equation for the evolution of M^f remains unchanged (neglecting the diffusion of Li 
through the sides of the cylinder for simplicity) . 

Unfortunately, the addition of a diffusive component to the Li flux now prevents us 
from deriving simple analytical solutions of these equations analogous to the ones found 
in §5.2.21 Instead, solutions can only be obtained numerically by integrating the partial 
differential equation fH2|) with time, and following the evolution of the spatially varying Li 



^Note that given all other approximations made in this section, the effects of radiative levitation and 
gravitational settling are neglected for simplicity. 
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profile. Given the other approximations made throughout this work (e.g. assuming that the 
rotation rate of the star, the background stellar structure and the convective forcing are all 
constant with time), and since our aim is merely a preliminary "proof-of-concept" , it seems 
futile to try to obtain a precise numerical solution of 



Instead, in this first paper we proceed by approximating the gradient term to cast (H2l) 
in the form of an ordinary differential equation similar to (!39|) . This can be done simply be 
writing 

VXu = ^^^:^"-° , (44) 

"diff 

where we have assumed that there is a "reservoir" of material with primordial Li abundance 
in the radiative zone in the form of a spherical shell of width d^is adjacent to rout, and 
that Li has to diffuse across that reservoir to be mixed into the convection zone. For this 
reservoir to contain enough Li to replenish the convection zone it must have roughly the 
same mass, so we set d^is = (3dout where /3 is a free parameter of the problem of order unity. 
Finally, the diffusion coefficient must also be approximated by a constant for (1421) to become 
a true ODE. Since the turbulent transport is weakest furthest from the radiative-convective 
interface, this is the "bottleneck" region for the diffusion of Li back into the convection zone. 
Hence we take 

D = -Dmicro(?^out " C^diff) + V°^^^dout exp ( -"T^ I . (45) 



A, 



As before, we divide 1^ by M°^^ to get 



where 
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^Lf - V^ - 2^ir^ ' (46) 

Li, pump Li, pump 

^ D 47rrg^tC^outp(r-out) _ D , .„^ 

"''' rfdiffrfout M°f /3t/L ■ 

As shown in Appendix B, one can integrate these equations analytically fairly straightfor- 
wardly. We find, as expected on physical grounds, that when the Li diffusion timescale into 
the convection zone becomes shorter than the advection timescale out of the convection zone, 
the surface Li abundance remains close to the primordial value. 

The resulting depletion fractions calculated using this method, for Li and Be, are shown 
in Fig. [15] for a simple grid of parameters (02 = 0.005, 02 = 0.02, (3 = 0.9 and (3 = 1.1), 
with the stellar models otherwise exactly the same as the ones used to create Figure 
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Note how, for all chosen parameter values, the Li and Be profiles now exhibit a clear dip 
centered roughly around the position of the observed Li-dip. On the cool side the effects 
of diffusion are negligible and the predicted depletion rates are very similar to the ones 
obtained in the absence of diffusion (see Figure [T^ . On the hot side, diffusion is important 
and continuously replenishes the outer convection zone with Li-rich and Be-rich material. 
The transition between the cool and hot sides, in the model, occurs when the diffusion 
timescale of Li back into the convection zone, namely Tdifj, becomes comparable with the 
advection timescale of Li-rich material out of it (namely T£?p^j^p), and depends both on the 
diffusion rate (as controlled by /3 for example) and on the pumping rate (as controlled by 
as). 

Note that for the larger value of 02, the ratio of the Li depletion fraction to the Be 
depletion fraction within the dip remains close to 1, contrary to observations. Instead, the 
ratio of the predicted depletion fractions is closer to observations for lower values of 02. The 
overall "best fit" is found for 02 = 0.005 and (3 = 0.9, which are not unreasonable parameter 
values. However, with our very crude model of the diffusion, we are not able to fit the exact 
position and amplitude of both Li and Be dips simultaneously. Since we were only aiming 
for a proof-of-concept, we view the results obtained as quite satisfactory, although a more 
careful model will be required in the future should one wish to explain the structure of both 
dips in more detail. 



5.4. Summary 

To summarize this section, we have found that the combination of gyroscopic pumping 
and turbulent diffusion of chemical species by overshooting motions ubiquitously predicts the 
presence of a dip in both Li and Be surface abundances for young MS stars in the mass range 
1.3 — I.5M0 (Li-dip stars). Depletion profiles close to the observed ones can be reproduced 
for reasonable values of model parameters. The increase in the depletion fraction on the cool 
side of the dip is explained by gyroscopic pumping and the progressively smaller amount of 
material which needs to be recirculated, while the decrease in the depletion fraction on the 
hot side of the dip is explained by the effect of diffusion on replenishing the outer convection 
zone with Li and Be. 
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Fig. 15. — Relative depletion of Li (solid lines) and Be (dashed lines) in the Hyades cluster 
as a function of stellar mass in the Li dip, as predicted for various values of the model 
parameters 02 (related to the stellar differential rotation) and (3 (related to the diffusion of 
chemical species across the base of the convection zone). AA^ is defined as in Fig. [131 See 
main text for detail. 



6. Discussion and prospects 



In this paper we performed an exhaustive study of a non-local source of rotational 
mixing called "gyroscopic pumping", which was originally studied in the context of the 
Earth's atmosphere by Haynes et al. (1991) and later discussed in the case of the Sun by 
Gough & Mclntyre (1998), Mclntyre (2007), and GAA09. 

In this mechanism, large-scale meridional motions are driven by angular-momentum 
conservation whenever fluid undergoes forces in the azimuthal direction (see §1.21 and in 
particular Figure [1] for detail). This is exactly the case in stellar convective zones, where 
rotationally influenced convection gives rise to the so-called A— effect (cf. Riidiger 1989; see 
also Garaud et al. 2010), and typically leads to the azimuthal acceleration of equatorial 
regions and deceleration of the poles. The resulting large-scale meridional circulation has 
fluid flowing outward from the rotation axis near the equator, and toward the rotation axis 
near the pole. Note that this mechanism is quite generic, since it is simply based on angular 
momentum conservation. However, by contrast with other well-studied signiflcant sources of 
rotational mixing, it does not rely on stellar spin-down to be effective - it is an inherently 
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quasi-steady mechanism. 

As first shown by GAA09 and studied more extensively here, whether the fluid gyro- 
scopically pumped in the convective zone ends up mixing the nearby radiative zone or not 
depends on many factors. The case of stars with a single convection zone, in the absence 
of any other mechanism, was first studied by GAA09. In this case the overall amplitude 
of the flows penetrating into the radiative zone is limited to slow Ekman flows, which are 
unlikely to play any significant role in mixing the stellar interior. However, this conclusion 
could be different if a large-scale magnetic field is present (Gough & Mclntyre 1998), as it is 
for instance thought to be the case in the Sun. We defer the magnetic case to a subsequent 
paper. 

Here we focused on stars with two convection zones, in the absence of magnetic fields 
(we accept that this approximation is probably over-simplistic). We found that gyroscopic 
pumping provides a significant source of non-local mixing in the radiative zones of these 
stars, more precisely along the rotation axis, within the cylinder tangent to the inner core 
(see Figure [TT]). This mixing readily explains, for example, the cool side of the well-known 
Li-dip and, when moderated by the effects of diffusion, provides predictions for the Li- and 
Be-depletion fractions which give a surprisingly good fit to the data on both sides of the dip 
(see Figure [T5|) . It is important to note that contrary to previous models, turbulent mixing 
here is needed not to explain Li destruction but to explain the replenishment of the surface 
layers in Li. Of course, the various other effects described in §5.1^ which were not taken into 
account here, could also affect the Li and Be abundances: mixing by gravity waves, large- 
scale rotational mixing induced by stellar spin-down, radiative acceleration, etc. A more 
sophisticated study of the Li-dip in the light of gyroscopic pumping, taking these effects into 
account as well as stellar evolution, is deferred to a subsequent pubhcation. Furthermore, it 
is interesting to note that the gyroscopic pumping mechanism probably plays a role in the 
global redistribution of angular momentum within the star. Whether this is related to the 
dichotomy between slow rotators and fast rotators on either side of the dip is an interesting 
question, which we hope to address in the future. 

Gyroscopic pumping, as studied here, also opens up other interesting observational 
prospects beyond Li-dip stars. For example, another important evolutionary phase when 
stars have two convective zones is the core helium burning phase on the red giant branch. Our 
study shows that this pumping mechanism could provide an important connection between 
the convective envelope and the core, and perhaps help reconciliate some of the long-standing 
discrepancies in surface abundances of CNO elements between models and observations. It 
is also interesting to note that the pumping itself has a tendency to drive a strong azimuthal 
shear in the system, close to the tangent cylinder (see Figure [6]). This shear could perhaps 
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become unstable for very rapidly rotating stars, leading to much more violent mixing in the 
radiative zone. Alternatively, this shear layer could also become the seat of an fi— effect and 
generate strong toroidal fields around the tangent cyhnder. It is interesting to speculate on 
what kind of field such a flow structure would sustain, and what observable feature it may 
lead to - although a full answer to that question will require a study of gyroscopic pumping 
in the presence of magnetic fields. 
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Appendix A: Analytical derivation of W^z 

In this appendix, we derive the analytical expression for the solution of the Cartesian 
model described in §3.2. ![ focusing in particular on the derivation of the mass flux into the 
radiative zone, Wj-z- Solutions to the set of equations ( 1T^ are first found in each of the three 
regions separately, then matched to one another at the interfaces (at ^in and ^out) and to 
boundary conditions (at z = and z = 1). 

Solution in the outer convection zone. In the outer convection zone, the viscous 
stresses are negligible compared with the linear drag term and N{z) =0. As a result, the 
equations reduce to 

2u = -ike'^^^pp - Aoutt' , 
= -p.e^/^" + e^/^^t - Aoutt^ , 

e-'/^'ikv + [e-'l^'w), = . (48) 

The temperature equation is easily solved as 

f{z) = ae^"" + he-^' . (49) 
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The remaining system of four equations can by transformed into a linear system of ODEs 
with constant coefficients using the following new mass flux variables: U = e~^^^'"u (and 
similarly for Ucz), V = e~^^^pv, and W = e~^f^pw. In that case, 

-2v = -Aout(f/-t/r), 

2U = -ikp - K^tV , 

ikV + W, = Q , (50) 

where L"^ = D'^ - D^\ 

Eliminating each variable in turn we can reduce the system to a second-order, forced 
linear ordinary differential equation for W (for example): 

W.. = W^^ + f^^r - e-^/^/^T , (51) 

4 + AL 4 + A2^t 4 + A2„t 

where 

STiz) = ^^^ = ^ (e-^/^^nr(z)) . (52) 

The solution of this equation is the sum of the solution of the homogeneous problem plus a 
particular solution. If we assume for simplicity that 5*°^* is constant in the outer convection 
zonci , the particular solution is easily expressed as: 



9,-C'out Ak-L i)z t J-k-L ^)z 

W{z) = Ae^l'-^+Be-"-^ ^'^^ a e be 



'cz 



Mout Aout 5l,,{k - L-i)2 - 1 Aout 5l,,{k + L-i)2 - 1 ' 

(53) 
where we have defined the new lengthscale 

. _ V4 + Ag„t .. .X 

°"* ~ A-A . ■ ^ ' 

Using the solutions obtained for W and T we deduce the solution for the pressure perturba- 
tion: 

p{z) = _^-Ao,t5out(^e^/'°"'-5e-^/'°"') 

{k - L'^)e^^-^''> (A; + L-i)e(-^-^"')^ 



^Of course, this is not the case for the example chosen here. However, the small difference in the solution 
caused by the non-constant 5™' is not worth the increased complication in the algebra since we are looking 
here only at obtaining a physical intuition of the solution. 



-45- 



Solution in the inner convection zone. By analogy, the solutions in the inner convec- 
tion zone are given by 



W{z) = Ce"/"^- + De""/^- - 



2^ 
kAw 



^{k-L-^)z 



d 



^(-k-L~^)z 



Ain S^^{k - L-i)2 - 1 Ai, 5l{k + L-i)2 - 1 



p{z) 



2im 



- Ain^m (Ce"/^'" - De-'/^'^) 



+c 



{k - L-^)e^^-^~"> Ak + L-^)e^~^-^'"> 



{k-L-^Y-5, 



-2 



{k + L~^] 



5 



-2 ' 



(56) 



where 5\-a and S™ are defined by analogy with 5out and 5°"* and where four new integration 
constants (c, rf, C and D) have been introduced. 



Solution in the radiative zone. In the radiative zone, the governing equations flT6|l 
reduce to 

-2v = Ey [u^^ - k^u) , 

2u = -ikpe'^^" + E^ {v^z - k'^v) , 

= -p.e^/^" + f e"/^^ + E, {w,, - k'^w) , 



-e ''-^'w = i^^\i^^-k T I , 



ikve-'^^' + {e-'l^'w), = . (57) 

In the limit of very small Ekman number, and for a < 1, the system is approximated by 

-2v = 0{E^ 
2u = -ike'/^^p + 0{E^ 






0{E, 



fe-^/^-u') 



(5^ 



which successively implies that 



^kz 



kz 



T = Kie"' + K2e-''' + 






W,, + 0{E, 



oz/L 



^l E, L-2 

2 



A;2 



P 



X (,z{k-L 1) j^ z{~k-L 1) JY2 ly 
Prz + -^ — . . . . + 



k-L-^ 



k + L 



-1 



ni E^L 



-2 



A;2 



+ 0{E, 



+ 0{E, 



(59) 
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and where IVrz, Ki, K2 and prz are integration constants. 

We see that the condition where pw is constant arises simply from mass conservation 
and geostrophy. Assuming that Wj-^ is known, the pumped flow induces local temperature 
perturbations as a result of the thermal energy equation: this can easily be seen in the 
solution for T, which contains two parts: the particular solution of the Poisson problem 
which includes the source term arising from the advection of the background temperature 
by the meridional flows, plus the general solution of V^T = 0. The constants Ki and K2 are 
simply there to match the temperature profile in the radiative zone to that of the convection 
zones. The expression for the pressure perturbation can be interpreted in a similar way. 

Matching the solutions and radial velocity in the radiative zone. We have found 
twelve integration constants, which can be uniquely determined by applying four boundary 
conditions (impermeability a.t z = and z = 1, and T = at the top and the bottom) and 
eight matching conditions (continuity of w - alternatively oi W - p, T and dT/dz at each 
of the two interfaces) . 

The boundary conditions aX z = and z = 1 imply: 



Mout Aout 6l^,{k - L-i)2 - 1 A„,t 5l^,{k + L-i)2 - 1 

r. r. 2^5"^ C 1 d I 

C + D = -r-^ + 



Mi, Ai^6l{k-L-^y-l Ai^6l{k + L-^y-l ' 

ae^ + be''' = , 

c + d = 0. (60) 

Continuity of vertical velocity, pressure, temperature and derivative of temperature at 
the z = Zont interface implies 

/l„2out/<5out _|_ T^ -Zout/Sout _ '-^cz [f '^ [^ '^ _ yrr 

A;Aout Aout 6l,ik - L-^Y - 1 Aout 5t,{k + L-^f - 1 

= Prz + 



k - L-1 k + L-1 f]2 ^^ 1^2 __ J.2 



JU2 TI/ pZout/L 

Qgfc^out ^ ^g-fczout ^ Kie'''°^' + 7^26"^^°"' + '' "' 



ni E^ L-2 - A;2 ' 



/V2 W pZout/L 

A;ae'=^-' - /tfoe"^"-' = kKie'''°-' - /tiTse"^"-' + -^— ^-- . (61) 
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Similar matching conditions at the lower interface {z 



imply 






2iSu 
kkw 



,{k^L-^)z;^ 



,{-k~L-^)z;^ 



Ai„ 5l{k - L-i)2 - 1 Ai, 5l{k + L-i)2 - 1 



W, 



rz ) 

-1^ 



k 



- 5inAin [Ce''-'^'- - De-"'"/'^'"] + 






s; 



d- 



{k + L- 



n2 



6Z 



Prz + ^ — ^-^-, + 



k-L- 



k + L- 



ni E, L-2 - fc2 ' 









1]2 R L-2 - A;2 ' 



o^^i 



/ti^ie'^"- - kK2e~'''''' + 



kZjYi 



Nl Wr: 



o^in/L 



n^ LE^ L-2 - A;2 



(62) 



This system of twelve equations can in principle be solved analytically exactly, but the 
solutions are horrendously complicated and without much interest. To get a better physical 
intuition of the solutions, we restrict our study to stars where the convection zone depths 
din and rfout are small compared with k, 6i^ and 6out- Since k = 2, and 5in and 6out are always 
greater than one for all possible value of Ajn and Aout, it is sufficient to require that rfin and 
dont be much smaller than one. Note that the din ^ 1 and dout ^ 1 approximations are 
acceptable for all Li-dip stars, which have thin outer convective regions and small convective 
cores. We also note that L is typically small compared with 6 or 1/k, although could be of 
the same order of magnitude as d-^i and dout- 

We begin by solving for the constants a, b, c, d, Ki and K2 in terms of W^z- This yields, 
in the limit of thin convection zones: 

'^ 2ni E^ V^kL 



Nl W„L^ 



b = K2+ 'I " e^out/Lgfc^out I I 



2fi2 R 



1 



kL 



,^ Nl W,zL' ,. n _k,. f, 1 

c = Ki-\ ^^ — e'^'^/^e ^"^'^ [l-\ 

'^ '^ kL 



d = Ko 



29.1 E^ 
NlW^ 

2fi2 R 



pZiii/ L/ kziii 



and then 








K, -- 


NlW,z L^ 


gZout/Lgfe f ij^ "out \ 

V L J 


ni E,i- 62^^ 


K2 -- 


Nl W,z L^ 






ni E^l- e-2*^ 



kL 



oZin/L 



(63) 



1j 



- e^'"/^ 1 - 



Win 

T 



(64) 
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Next, we solve for A, B, C and D in terms of VFrz: 

"^^out ^out Oout V "out/ ^ "out ^outOout"out 



5g-i/5o.. ^ ^^ _ ^-z^u^/L^kl^ ^ _ ^^ ^ M;^^ _ ^ae^fc^^_i/L 



"^Aout Aout Oout V "out/ 2 dout AoutOout"out 



^ — i~i T ^ ^^ "; ;^ — I -L H — ; — ~r 



rvlVirt -^i-in '-'in \ (^\r\ / ^ Z\n Z\x^l\.\r\0\ 



D = ^- e-WL^^ f 1 + - V — - + ^'^^[ ' (65) 



where a and c are explicitly written above. 

These expressions can finally be substituted into the pressure continuity equations at 
both interfaces, yielding two equations for p^z and W^z- Eliminating pr^ yields (JTOll . A gener- 
alization of this derivation for more complex density, temperature and buoyancy frequency 
profiles can also be performed analytically although in practice is not particularly useful. For 
real stellar models, it is much easier to solve (ITB]) numerically, although the overall scaling 
of the solutions with parameters remains the same. 



Appendix B: Li and Be mass fractions 

In this appendix we solve the evolution equation for the Li abundance in the outer 
convection zone with and without diffusion. The same calculation can be done to evaluate 
the evolution of the Be abundance. 

The set of equations we wish to solve, (1391) or (H6|) . can be cast in the following general 
form: 

Xcz v'^y^ VC7, V 
Alj — 1 , 

rj-cz ^cz T^-fr 

pump pump oirt 

^Lf = -^^ - ^^^ ■ (66) 



cyl cyl 

Li, pump Li, pump 

To reproduce equation ( 15^ we set rdig -^ oo. Initial conditions are 

^Li(0)=X^f(0) = XLi,o. (67) 



This set of equations has the solution 



A^Li.O 



Tdiff 
pump 



+ 1 ' 



-49- 



-cyl/ 



X^'it) = 1 + AiC„,p + 



r 



pump 
Tdiflf 



Ae^^* + 1 + Aar; 



r. 



pump 



pump 
Tdiff 



5e 



A2t 



X, 



Li,0 



Tdiff 



^ + 



■63) 



where 

Al,2 = 



^cyl 



Li,pump pump ""i 



± 



\40 



-cz A 2 
pump / 



^''~pump''~diff 



^cyl 



^Li.pump 2rdiff 



(69) 



Applying the initial conditions yields 



A = Xu, 



A, 



A2 — Ai 1 -|- 2 ^p"'"p 

Tdiff 



-B = Xli n 



Ai 



Ai — A2 1 + 2 ^p""'p 

Tdiff 



(70) 



Let us first examine the case with no diffusion (rdig — ?> 00). In this case 



y{t) = (1 + Air-^p) Ae^^* + (l + X^r^l^^) Be'-' 



with 



Ai,: 



1 



1 



^cyl 



2^cz 
pump 



± 



4(r, 



cz ^2 

pump/ 



+ 



1 



^cyl 



(71) 



(72) 



'Li.pump -'pumpy y ^Vpumpy V'Li.pump^ 

In the limit where the mass of the outer convection zone becomes small compared with 
the mass of the cylinder Cu-, then Tp^^^p — > 0. This implies that X^^{t) ~ X^^ (t), so 
that the respective concentrations of Li (or Be) in the convection zone and in the cylinder 
are the same. Moreover, they vary essentially as e~*''^Li,pump since when Tp^^^^p <ti T^^nnrapi 
Ai ~ "l/'^'upump ^"^^ A2 — — 1/Tpump- The evolution on the timescale Tp^^^p is extremely 
fast, leaving the longer timescale TLfp^jj^p as the depletion timescale for the system. In the 
opposite limit (with the mass of the cylinder going to 0) it is easy to show that as before 
X£;(t) ~ X^f (t) and the depletion timescale is now Tp^^^^p. 

We expect that the effects of diffusion become important when the diffusion timescale 
into the convection zone Tdiff becomes shorter than the advection timescale out of the con- 
vection zone Tpyjjj . In the limit Tdiff/Tp^ ^ 1 then the solution approximately becomes 



Xli.O 



Tdiff -t/r"?'' 

p ' iji,pump 



rl- 



diff 



x'^m 



2t-cz 
pump 



p ' Li, pump 



o^cz ^cyl 

pump Li, pump 



-t/Tdiff 



'diff 



^cz ^cyl 
pump Li, pump 



-*/Tdiff 



+ 1 



XLi,0 



XLi,0 ; 
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so that the Li abundance in the convection zone remains roughly primordial, as expected. 
The Li concentration in the cylinder gradually tends to half the primordial abundance, which 
is also expected since the cylinder itself becomes fully ventilated, and is a mixed region with 
equal incoming mass flux with primordial Li abundance, and incoming mass flux with zero 
Li abundance. 
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M, 


R. 


din 


"out 


conv 


„.out 
conv 


ru 


^Bc 


fi. 


1.30 


8.84 


5.5 X 10-2 


1.2 X 10-1 


8.0x102 


2x 10^ 


0.61 


0.51 


8.8 X 10-5 


1.35 


9.28 


5.8 X 10-2 


1.0 X 10-1 


5.0x102 


1 X 10^ 


0.60 


0.50 


1.0 X 10-4 


1.40 


9.69 


6.3 X 10-2 


7.5 X 10-2 


1.0 xlO^ 


3x 10^ 


0.59 


0.49 


1.3 X 10-4 


1.45 


10.1 


6.9 X 10-2 


5.3 X 10-2 


1.2 xlO^ 


5x 10^ 


0.57 


0.48 


1.5 X 10-4 


1.50 


10.3 


7.5 X 10-2 


2.9 X 10-2 


1.5 xlO^ 


1 X 10^ 


0.57 


0.48 


1.5 X 10-4 


1.55 


10.6 


8.0 X 10-2 


1.4 X 10-2 


1.8 xlO^ 


1 X 10^ 


0.57 


0.48 


1.5 X 10-4 


1.60 


10.8 


8.5 X 10-2 


6.5 X 10-3 


2.0 xlO^ 


1 X 10^ 


0.57 


0.48 


1.5 X 10-4 



Table 1: Stellar features as extracted from our reference 1.3M0 — I.6M0 background models 
evolved to 300 Myr. From left to right: the stellar mass (in units of Mq), the stellar radius 
(in units of lO^^cm), the depths of the inner and outer convection zones in units of the stellar 
radius, the convective velocities in each convection zone in cm/s, the Li and Be burning radii 
(see main text) in units of the stellar radius, and finally an estimated rotation rate for the 
star (see main text) in radians per second. 



